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Abstract 

In this work, we discuss the implications of a recently obtained equi- 
librium fluctuation-dissipation relation on the extension of the available 
Monte Carlo methods based on the consideration of the Gibbs canonical 
ensemble to account for the existence of an anomalous regime with nega- 
tive heat capacities C < 0. The resulting framework appears as a suitable 
generalization of the methodology associated with the so-called dynami- 
cal ensemble, which is applied to the extension of two well-known Monte 
Carlo methods: the Metropolis importance sample and the Swendsen- 
Wang clusters algorithm. These Monte Carlo algorithms are employed to 
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study the anomalous thermodynamic behavior of the Potts models with 
many spin states q defined on a d-dimensional hypercubic lattice with 
periodic boundary conditions, which successfully reduce the exponential 
divergence of decorrelation time r with the increase of the system size 
to a weak power-law divergence r oc A''" with a ~ 0.2 for the particular 
case of the 2D 10-state Potts model. 

1 Introduction 

In the present work, we shall not propose new Monte Carlo (MC) methods based 
on the equilibrium distributions of Statistical Mechanics. On the contrary, we 
shall discuss how the available MC methods based on the consideration of the 
Gibbs canonical ensemble: 

dp {E \I3b ) = Z {PbY^ exp {-PbE) n (E) dE (1) 

could be extended by using a minimal, but crucial modification in their schemes 
to account for the existence of an anomalous regime with negative heat capac- 
ities C < [TJ[U[31[11[5J[B]. This fact avoids the incidence of the so-called 
super- critical slowing down, a dynamical anomaly that significantly affects the 
efficiency of large-scale canonical MC simulations [7] . 

Our proposal follows as a direct application of the recently obtained fluctuation- 
dissipation relation [SI H] : 

C = I3^{SE'') + C{6ME), (2) 

which involves the heat capacity C of a given system in an equilibrium situation 
where the inverse temperature /3(j of a certain environment exhibits correlated 
fluctuations with the system internal energy S as a consequence of their mutual 
thermodynamic interactiorQ. Eq.© accounts for the realistic possibility that 
the internal state of the system acting as environment could be affected by 
the presence of the system under study, which is a fact a priory disregarded 
by the consideration of the Gibbs canonical ensemble ((T|), where the inverse 
temperature /3i^ of the environment exhibits a constant value fis because its 
heat capacity Ci^ is practically infinite. Obviously, Eq.® is just a suitable 
extension of the well-known relation: 

C^P^ (SE^) (3) 

between the heat capacity C and the energy fluctuations derived from the canon- 
ical ensemble H]). While the canonical result ^ only admits macrostates with 
positive heat capacities C > 0, it is easy to verify that the fluctuation relation 
^ is compatible with the presence of macrostates having negative heat capac- 
ities C < 0. This last conclusion is the fundamental ingredient considered in 
this work for allowing a direct extension of some MC algorithms based on the 

^ Along this work, Boltzmann constant is assumed to be fcg = 1. 
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Gibbs canonical ensemble (O in order to account for a regime with C < 0. As 
discussed elsewhere [U [2 [3l IH [5l [6] , this kind of anomaly in the caloric curve 
appears to be associated with the occurrence of a discontinuous (first-order) 
phase transition in finite short-range interacting systems, as well as systems 
with long-range interactions such as astrophysical systems. 

This work is organized into sections as follows: first, we shall discuss in sec- 
tion [2] how the present ideas can be considered to extend the available canonical 
MC methods; afterwards, we shall apply these arguments in section [3] to extend 
two well-known canonical MC algorithms: the Metropolis importance sample 
|101 111] and the Swendsen- Wang cluster algorithm [121 1131 114] , in their appli- 
cation to the study of anomalous macrostates present in the thermodynamic 
description of the g-states Potts models defined on a d-dimensional hypercubic 
lattice; finally, concluding remarks are presented in section |4l 



2 The proposal 
2.1 Overview 

For convenience, let us begin the present discussion by reviewing the most im- 
portant results related to the fluctuation theorem ([2]). Our analysis starts from 
the consideration of the following generic energy distribution function ^ : 

dp^ (E) ^uj{E)n (E) dE, (4) 

where w (E) is a probabilistic weight that considers the thermodynamic influence 
of a certain environment. The above work hypothesis admits the canonical 
weight: 

iJciE) = Z{/3By'ex_p{~l3BE) (5) 

as a relevant but particular case when the environment is just a thermal bath 
having an inflnite heat capacity. This general situation could be implemented 
with the help of a Metropolis Monte Carlo simulation by using the transition 
probability: 

WAE^E + 5E) = min | "(f + ^^\ l} . (6) 

Since the energy thermal fluctuations \SE\ are small when the system size is 
sufficiently large, Eq.([6]) can be rewritten in a canonical fashion as follows: 

W^{E^E + SE) ~ min {cxp [-/^^ (E) SE] , 1} , (7) 

where (E) is hereafter referred to as the inverse temperature of the environ- 
ment: 

/3.(i^)^^ = -4log.(i5). (8) 

The density of states ft (E) is related to the system entropy S (E) as ft (E) = 
C exp [S {E)], which allows us to obtain the system temperature T by using the 
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thermodynamic relation; 



1 _dS{E)_ 

Eqs.(|8l) and ([9]) can be combined to express the inverse temperature difference 
Tj as follows: 

il{E)^P^{E)^P{E) = -^\ogp{E), (10) 

with p (E) — uj (E) n (E) being the density of probability. This last representa- 
tion (|10p is very useful to obtain two remarkable thermodynamic relations. The 
first one involves the statistical expectation value (rj) , and its calculation reads 
as follows: 

{ri)=j^ v{E)p{E)dE = - Q^p{E)dE 

= -p(^)lt:r =0, (11) 
while the second one considers the correlation function {Erj) : 

{Ev) - ^ Erj (E) p{E)dE^-J^ E-^p {E) dE 

p{E)dE = l. (12) 

-Eint 

Here, we have taken into account the vanishing of the density of probability 
p (E) and its first derivative dp (E) / dE at the maximum Egup and minimum 
ii'inf values of the system energy, as well as the normalization condition. 

The vanishing of the expectation value (ry), Eq. (jll[) . is simply the known 
thermal equilibrium condition: 

the mathematical form of which clearly indicates that the equalization of tem- 
perature expressed by the Zeroth Principle of Thermodynamics takes place in 
an average sense. Eq.()12|) can be rewritten as a rigorous fluctuation relation by 
using the identity {dEdrj) = (Erj) — {E) {rf) : 

By using the Schwartz inequality (AB)'^ < (A^) (i?^), this last result can be 
rephrased as: 

AEA - ^) > 1, (15) 
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where Ax — ^ {dx^) . Finally, by substituting the first-order approximation: 



into Eg. ([T4| . with C = dE/dT being the heat capacity, one obtains the fluctuation- 
dissipation relation ©. 

The fluctuation-dissipation relation ([2]) can be rewritten as follows: 

C{l-{6(3^5E))^p^{5E''). (17) 

Since the right-hand side of this expression is always nonnegative, it is easy 
to see that the presence of macrostates with positive heat capacities C > 
demands that the correlation function {5P^5E) obey the constraint: 

{5PJE) < 1. (18) 

Clearly, such a condition is fulfilled by the Gibbs canonical ensemble ([T]), where 
5(3^ = 0. However, the existence of macrostates having negative heat capacities 
C < can be only observed provided that the constraint: 

{SME) > 1 (19) 

holds. Thus, any attempt to impose the canonical condition (5/3(^ — ;> is always 
accompanied with a progressive increase of the energy fluctuations JE' — >■ oo, 
which leads to the thermodynamic instability and inaccessibility of such anoma- 
lous macrostates. 

The simplest way to guarantee the existence of non-vanishing correlated 
fluctuations {SPi^SE) 7^ is achieved by considering an environment with a 
finite heat capacity C^o- Here, the inverse temperature fluctuations SP^j can be 
expressed in terms of the amount of energy 5E released or absorbed by the 
system in turn its equilibrium value: 

Sp^^p-'^SE, (20) 

where we have considered the thermal equilibrium condition /3[^ = By substi- 
tuting this last expression into fluctuation-dissipation relation ([2]), one obtains: 

^^'^ =P^5E^). (21) 



Since the right-hand side of this last expression is always nonnegative, the ther- 
modynamic stability of macrostates with negative heat capacity C < demands 
the applicability of the following constraint [9] : 

Q<C^< \C\. (22) 

Remarkably, this result was also obtained in the past by Thirring |15) . 
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Figure 1: Schematic behavior of the microcanonical caloric /3 (E) — dS {E) /dE 
of a finite short-range interacting system undergoing a first-order phase transi- 
tion between a low (lep) and a high (hep) energy phase. The density of proba- 
bilities pi (E) and p2 (E) are associated with the thermal contact of this system 
with certain environments characterized by the inverse temperatures (3^ (E) and 
(3'^ (E) respectively. 



The above consequences are illustrated in detail in FIGlH We show here 
the typical backbending behavior of the microcanonical caloric curve (3{E) of 
a finite short-range interacting system undergoing a first-order phase transition 
[3] , where the points within the energetic region Ep < E < Eq represent anoma- 
lous macrostates with negative heat capacities. The density of probability p (E) 
corresponds to a situation where this system is put in thermal contact with a cer- 
tain environment characterized by the inverse temperature (E) . The maxima 
and minima of the distribution function p (E) are determined from the thermal 
equilibrium condition /3 [E) — 13^^ (E), that is, the intersection points between 
these inverse temperature dependencies. For convenience, we have shown here 
two relevant cases. 

The first case corresponds to the equilibrium situation associated with the 
Gibbs canonical ensemble ([TJ , where the inverse temperature dependence /3]j (E) 
remains at the constant value /3b despite the underlying energy interchange. It 
should be noticed that the thermal equilibrium condition is fulfilled by three 
points {Ea, Eh, Ec) when the inverse temperature /3b is within the interval 
{l3p,l3q). If this is the case, the energy distribution function pi (E) is bimodal, 
where the points Ea and Ec determine the positions of its peaks (local maxima), 
while the intermediate point Eh determines the position of its local minimum. 
One can verify that the local minimum Ef, always belongs to the anomalous re- 
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gion with C < 0, while the maxima {Ea,Ec) are located within the regions with 
C > 0. The accessibility of such points behaves with the increase of the system 
size N as pi [Ea^c] oc e"°'=^ and pi (Et) oc e""'^, where aa,b,c > 0. Con- 
sequently, anomalous macrostates with C < becomes practically inaccessible 
within the canonical ensemble ([1} when iV is sufficiently large. The existence of 
such a hidden region is the origin of the latent heat ql necessary for the conver- 
sion of one phase into the other during the phase transition, as well as for the 
ensemble inequivalence between the microcanonical and canonical description 



A multi-modal character of the density of probability pi (E) in the framework 
of MC simulations based on the canonical ensemble ([1]) leads to the occurrence 
of the super- critical slowing down [7]: an exponential divergence of correlation 
times with the increasing of the system size, r oc exp (AiV). In other words, this 
phenomenon manifests as an effective trapping of the system energy within any 
of the coexisting peaks of the distribution function pi {E), due to the probability 
T for the occurrence of a large energy fluctuation that allows a transition towards 
a neighboring peak decreases exponentially with the increase of the system size 
N, T ~ pi {Eb) I p\ (Ea) oc e^'*'^. Thus, the characteristic timescale for such 
a transition is given by r oc 1/T ^ exp (XN). Because the canonical averages 
should consider the contribution of all these coexisting peaks, the relaxation 
timescales of such expectation values also exhibit an exponential growth with 
the increasing of N. 

The second case shown in FIG H] corresponds to a situation where the system 
is put in thermal contact with an environment having a finite heat capacity. By 
choosing appropriately the environment and its internal conditions, in partic- 
ular, the applicability of Thirring's constraint (|22p. the corresponding inverse 
temperature f]^ (E) can ensure the existence of only one intersection point with 
the microcanonical caloric curve /3 (E) of the system under study. Even, such a 
point could be located within the anomalous region with C < 0, e.g., the unsta- 
ble macrostate _Ef,. Since the energy distribution function p2 {E) is monomodal, 
the phenomenon of super-critical slowing down cannot be present when such a 
physical situation is simulated by using a suitable MC method. 

2.2 Application in Monte Carlo methods 

In the multicanonical MC method and its variants , the main aim is to obtain 
of the density of states ft{E), or, equivalently, the microcanonical entropy S{E). 
Such a goal can be achieved through a direct MC calculation of the energy 
distribution function p{E) ~ lj{E)^1{E), which can be inverted to express the 
system entropy and the canonical expectation values (^)a as follows: 



S{E) = log [p{E)] - log [u^{E)] + cte, 



(23) 




uj{E) 



'^^dE 



dE 



(24) 
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The mathematical form of the probabihstic weight, uj{E), can be conveniently 
proposed before performing the MC simulation |16) , although the most common 
strategy is to carry out its iterative reconstruction through several preliminary 
MC runs to obtain a flat histogram p{E) ~ const within the energy interval of 
interest (TTHTB]. This latter alternative allows to enhance rare events, which is 
particularly useful to avoid the super-critical slowing down near to first-order 
phase transitions [7] . This kind of MC simulation allows the acquisition of the 
thermodynamic information in a wide energy interval by performing a single 
MC run. A clear disadvantage is that a substantial fraction of computational re- 
sources should be consumed to find an optimal probabilistic weight ijj{E). More- 
over, the acquisition of the microcanonical caloric curve (3{E) = dS{E)/dE as 
well as of the heat capacity C{E) — —(3^{E) [d'^S{E)/dE'^'\ from a numerical 
differentiation of the entropy S{E) is a procedure that enhances the unavoidable 
statistical errors involved in the MC calculation of the energy distribution p{E) 
and/or the probabilistic weight lu{E). 

Clearly, it is desirable to implement some type of methodology that allows 
a direct estimation of the caloric curve /?(£') and of the heat capacity C{E) 
without considering a numerical differentiation of the system entropy S{E). A 
simple strategy is provided by considering the thermal contact of the system 
with an environment characterized by a finite heat capacity, as in the second 
case illustrated in FIGHJ It should be noticed here that the system energy 
E and the inverse temperature of the environment undergo small thermal 
fluctuations around their expectation values (E) and which provide a 

suitable estimation of the intersection point {Eg, /3e) derived from the condition 
of thermal equilibrium Pui{Ee) = PiE^) = Pe- 

E, = {E) , /3e = . (25) 

In essence, Eq. (|25]) expresses the procedure associated with the MC method 
of the so-called dynamical ensemble proposed by Gerling and Hiiller [19j , which 
is defined by a probabilistic weight ijj{E) with a power-law shape: 

lude{E) = A^o {Et ~ Ef , (26) 

whose corresponding inverse temperature is given by: 

P^{E)^BI{Et~E). (27) 

The fluctuation relation ^ enables the introduction of several improvements 
for this kind of MC calculations to estimate the thermo-statistical properties of 
a system within the microcanonical ensemble by using modifled canonical MC 
algorithms. In fact, one can also obtain the heat capacity Cg — C{Ee) at the 
interception point from the fluctuating behavior of the system as follows: 

13^, (SE^) 

~ 1-{SME)- ^^^^ 

Thus, one is able to obtain a suitable estimation of any point on the micro- 
canonical caloric curve /3 (E) as well as the heat capacity C{E) regardless of 
whether its character positive or negative. 
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2.2.1 The linear ansatz and its optimization 

Generally speaking, the specific mathematical form of the probabilistic weight 
u}(E) or its corresponding inverse temperature Pcj{E) is unimportant as long as 
the following conditions applied: 

1. The inverse temperature Plo{E) of the environment must ensure the exis- 
tence of only one intersection point with the caloric curve /3{E), that is, 
the existence of only one peak in the energy distribution function p{E). 

2. The size N of the system under analysis should be sufficiently large to 
guarantee the validity of the first-order approximation (jl6l) considered for 
deriving the fiuctuation-dissipation relation ^ from the rigorous fiuctua- 
tion relation (ITil) . 

As already depicted in FIG HI the energy distribution function p{E) has the 
shape of a bell curve, usually approximated by a Gaussian distribution. Since 
the system exhibits small thermal fluctuations /S.E/\E\ ~ l/-\/iV, only local 
properties of the inverse temperature dependence Puj (E) close to the equilibrium 
point Ee are significant; hence, one can obtain the same practical results by using 
different mathematical dependencies for the inverse temperature l3u>{E). 

The simplest mathematical form is provided by restricting to the first-order 
power expansion of the inverse temperature dependence (E) around the in- 
tersection point: 

f3u{E)^l3, + X5E/N, (29) 

which considers a linear coupling of the environment inverse temperature /3(^ 
with the thermal fluctuations of the system energy SE — E — E^. Remarkably, 
such a linear ansatz (j29p is equivalent to the so-called Gaussian ensemble |20] : 



^Ge{E) = Aexp 



-/3ei? - ^KE - E,f 



(30) 



proposed by Hetherington [21 , which approaches in the limit A — !■ +oo to the 
microcanonical ensemble u!me{E) — AS{E — E^). 

Hypothesis (I^H]) ensures the existence of non-vanishing correlated fluctua- 
tions {5f5^5E) 7^ when the coupling constant A 7^ 0. By using the fluctuation- 
dissipation relation ([2]) and the ansatz (|29l) . one can obtain the following ex- 
pression for the heat capacity: 

C (31) 

which can be rewritten to obtain the energy and inverse temperature dispersions, 
AE and as follows: 

N 1 
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where Aa; = \/ {Sx^}. The inverse temperature dispersion A/3i^ decreases with 
the increase of the system size N as A/3i^ = \X\AE/N ~ 1/^/N. In the hmit 
A'^ — > +00, the thermal fluctuations disappear and the present equilibrium sit- 
uation becomes equivalent to the microcanonical ensemble provided that the 
following condition: 

l3lN/Ce + A > (33) 

holds. According to the first-order approximation (|20p . the coupling constant 
A can be related to the heat capacity C^^ of the environment as A = 0lN /C^^^ 
hence, the stability condition (|33|) is fully equivalent to Thirring's constraint 
(1221). 

At first glance, it is desirable to maximally reduce the thermal dispersions 
of the system energy E and its inverse temperature /3 = 1 /T indirectly derived 
from the environment inverse temperature /3(^ = 1/21,. However, the inequality 
of Ea.(|15p imposes an important limitation on the precision of such a kind of 
measuring process: a reduction of the thermal uncertainties A{\/T^ — l/T) 
affecting the temperature equalization provokes an increasing of the system 
energy fluctuations AE, and vice versa. As previously discussed [319], Eq. ([T5| 
accounts for the existence of some type of complementary character between 
thermodynamic quantities of energy and temperature. In consequence, one 
must assume the existence of non- vanishing thermal uncertainties for the system 
energy E and its inverse temperature /3. 

According to expressions (|32|) , the increase of the coupling constant A allows 
us to reduce the energy dispersion AE, but its value should not be excessively 
large because its increasing also leads to the increasing of the inverse temper- 
ature dispersion A(i^. A simple criterion to provide the best value for A is 
obtained after minimizing the total dispersion A'^ = {AEY / N + N {AjSi^Y . By 
introducing the microcanonical curvature k: 

K = k{E) = f3^N/C = -Nd^S (E) /dE^, (34) 

which is defined in terms of the second derivative of the entropy with the oppo- 
site sign, the optimal value (opt) for A is given by: 

Aopt (k.) = - K. (35) 

The microcanonical caloric curve /3 {E) and the heat capacity C (E) derived 
from Eqs. ([25l) and ([3T1) should not depend on the coupling constant A as long as 
this latter parameter fulfils condition (1551) and the system size N is sufficiently 
large. Nevertheless, the use of its optimal value ((35)) minimizes the underlying 
thermal fluctuations, which should also reduce the number of steps needed to 
ensure the convergence of a MC run. 

2.2.2 Iterative schemes 

Given a certain dependence of the environment inverse temperature /J^'^ (E) , one 
can obtain from a MC run suitable estimations of the system inverse temperature 
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j3i and its heat capacity Ci at the i-th interception point Ei. These estimates 
values can be employed to provide the next dependence (E) to perform an 

analogous estimation at the i + 1-st neighboring point To fix some ideas, 

let us denote by e = E/N the system energy per particle. The i-th dependence 
of the environment inverse temperature (pUj) considered for the MC calculation 
of the i-th point Si is given by: 

(e) = + A, (e - £-,) , (36) 

where and /3i are some rough estimates of the expectation values — (e) and 
Pi = ^/^i'^y The parameters {Si, Pi^Xi) could be provided in an interactive way 
by using the curvature Ki — /SfN/Ci derived from the previous MC simulation: 

Ei+i = £i + Ae, = A - KiAe, X,+i = Aopt (kj) , (37) 

where Ae is a small energy step. The initial value Sq is assumed to be the 
expectation value (e) of the energy per particle obtained from any canonical 
MC simulation with inverse temperature Pb = sufficiently away from the 
anomalous region with C < 0. 

The stability of the iterative procedure previously described depends on 
the precision of the estimated value of microcanonical curvature k. Of course, 
alternative iteration schemes are also possible, e.g., the following scheme: 

£-,+1 = (e)^ + Ae, /3,+i = , (38) 

which forces with Ae > (Ae < 0) a forward (backward) motion of the ex- 
pectation values (s) and (Pui) along the system caloric curve /? (e) whenever the 
value of coupling parameter A does not change in a significant way. In general, 
the inverse temperature /? (e) of a large enough short-range interacting system 
describes a plateau within the anomalous region with C < 0, which means that 
the corresponding curvature k does not significantly differ from zero, k ~ 0. 
In such cases, the value A ~ 1 constitutes a suitable approximation within the 
energy range containing the anomalous region with C < 0, while the value A = 
corresponding to the canonical ensemble is a good choice elsewhere. 



2.2.3 Implementation 

Since the inverse temperature Pb of the Gibbs canonical ensemble ([T]) appears 
as a driving parameter in the transition probability W{X — >■ X'; (3b) of canon- 
ical MC methods, the most general idea to extend this kind of algorithms is 
to replace the bath inverse temperature /3b by a variable inverse temperature 
Pui{E), Pb — ^ Pu{E). Moreover, it is desirable that the transition probability 
resulting from such a modification fulfills the so-called detailed balance condition 

p{X)W[X ^ X') ^p[X')W{X' ^ X), (39) 

where p{X) is the system distribution function, which is given by the function 

p^{X)^Lo[E{X)]. 
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Let us denote the energy varying during the configuration change X — > X' 
asSE = E {X') - E {X) and its mean value as E"" = [E (X) + E {X')] /2. Ac- 
cording to the mean value theorem, one can express the variation of the function 
logpui {X) by using the environment inverse temperature (3^1 {E) evaluated at a 
certain intermediate energy Eg = E"^ + 05E as follows: 

logp„ {X') - logp„ [X) = -P*JE, (40) 

where 9 is & real parameter with 16*1 < 1 and /3* = {Ee)- By considering the 
transition probability of any canonical MC algorithm that obeys the detailed 
balance condition: 

W{X^X';Pb) _ , . .... 
W(X-^X;/3,) =^"P^-^"^^)' (41) 

as well as Eg. piH) . one can finally obtain: 

Vu. (X) W{X^ X'; = {X') W {X' ^ X- /?*) . (42) 

This last result clarifies that given the initial and final system configurations, X 
and X', one can always find a certain value /3* = /J^j {X, X') of the bath inverse 
temperature that fulfils the detailed balance condition for the present dis- 
tribution function {X). In particular, the exact value of /3* for the Gaussian 
ensemble (15(1)) is the one corresponding to = 0, /3* = (3cj{E™'). 

The main obstacle to perform a direct application of Ea. (l42t to extend 
canonical MC methods relies on the fact that the final system configuration 
X' is a priori unknown in many non-local MC algorithms [HI [T31 HH HH [221 
[231 [23 [221 [23 . Consequently, one should consider some suitable approximation 
/3* = (Eg-) of the exact value /?*, e.g., the one corresponding to the environ- 
ment inverse temperature at the initial system configuration X, Pl^ = [i?(X)]. 
To estimate the error involved in this last approximation, it is important to take 
into account that the environment heat capacity C^o and the energy change 5E 
behave with the increase of the system size N as Cu: ^ N and 6E ~ N". Here, 
the exponent a ranges from zero for a local algorithm such as the Metropolis 
importance sample, up to 1/2 for a hypothetical non-local algorithm able to ob- 
tain an effective independent configuration after each MC stcfU. Consequently, 
the difference — — (3*\ merely constitutes a small size effect: 

S^-iPZ'f^m^j^, (43) 

where the index m indicates that the corresponding quantities l3i^ and C'l^ have 
been evaluated at the energy value E™ . Since the estimated inverse temperature 
depends on the system energy E(X), the corresponding distribution function 
associated with this approximation p^(X) can also be expressed by a certain 



^The exponent a = 1/2 follows from the size behavior of the energy dispersion AE ' 
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function of the system energy, p* (X) = uj* [E{X)]. According to Eg. ipS)) . the 
corresponding inverse temperature: 

,„.,^,._M) ,44, 

cannot differ in a significant way from the exact dependence f5i^{E) as long as 
the system size N be sufficiently large. 

The extended canonical MC methods based on an estimation of the inverse 
temperature P* do not fulfill the detailed balance condition ([39]). However, this 
fact does not represent any fundamental difficulty since one practically obtains 
the same numerical results for the caloric curve /3{E) and the heat capacity C{E) 
with the help of equations (|25|) and (|3T|) by using slightly different probabilistic 
weight uj*{E). The only requirement is that the energy distribution function 
p{E) exhibits a sharp Gaussian profile, which is simply achieved when the size 
N of the system under study is large enough. While the super-critical slowing 
down of canonical MC methods observed in systems undergoing a first-order 
phase transition becomes more severe as N increases, the errors associated with 
all approximations assumed here turn more and more negligible. This is the 
reason why the present methodology is particularly useful for avoiding this type 
of slow sampling problems in large-scale MC simulations. 

As naturally expected, finite size effects can be significant when one is also 
interested to describe systems whose size N is not so large. If this is the case, 
it is not only necessary to implement extended MC schemes that fulfills the 
detailed balance condition (|39|) . but also the inclusion of some finite size correc- 
tions into equations (|25|) and (|31|) employed here to obtain the microcanonical 
dependencies f3{E) and C{E). Although the complete analysis of these ques- 
tions is outsize of the scope of the present paper, we would like to clarify that 
the most general way to fulfill the detailed balance condition (p9| after the con- 
sideration of an estimated value for the transition inverse temperature /3* is to 
introduce a posteriori an acceptance probability w. 



w = mm 



l,^exp(-/3^5i?)^ (45) 



to accept or reject the final configuration X' . Here, Wi^f = W \X ^ X'\ 
and Wj^, = W{X' ^ with = p^{X) and = /3^(X'), represent 

the transition probabilities of the direct and the reverse process, respectively. 
The mathematical forms of which depend on the particularities of each non-local 
MC method. 
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3 Application examples 



3.1 The model 

For the sake of simplicity, let us consider in the present study the q-state Potts 
model [13]: 

defined on a d-dimensional hypercubic lattice N ^ Li x L2 x . . . Ld with Li = L 
and periodic boundary conditions. Here, the sum is only over pairs of nearest- 
neighbor lattice sites and cr^ = [1,2,... q] is the spin state at the i-th site. This 
model system is just a generalization of the known Ising model: 

^^Ising = - E SiSj , (47) 

with Si = ±1, which appears as a particular case with q = 2 after considering a 
linear transformation among their respective energy per particle £ising = — d + 
2£2 -Potts and inverse temperature /3ising = 5/32-Potts. 

As discussed elsewhere [7], this model system undergoes a first-order phase 
transition when g > 4 for d = 2 and g > 3 for d = 3, and hence, the mi- 
crocanonical caloric curves /3 (e) for these model realizations exhibit the back- 
bending behavior associated with the presence of macrostates with C < such 
as the one sketched in FIG[T] In addition to the consideration of a local MC 
method as the Metropolis importance sample \10\ lllj , the canonical MC study 
of Potts models can be carried out by using another accelerating methods such 
as the Swendsen-Wang [T^[T3] or Wolff [14 clusters algorithms. However, none 
of these MC algorithms are able to account for the existence of an anomalous 
regime with C < 0, and even, they also suffer from the existence of a super- 
critical slowing down as a direct consequence of the bimodality of the canonical 
energy distribution function. The existence of several canonical MC algorithms 
for this kind of models enables to perform a comparative study among their 
respective extended versions described below. 

3.2 Monte Carlo methods 
3.2.1 Metropolis importance sample 

The simplest and most general way to implement a thermal coupling with a 
bath at constant inverse temperature f^s is by using the Metropolis importance 
sample [inilll]. In this method, a Metropolis move is carried out as follows: 

1. A site i is chosen at random, and the initial spin state ct^ is changed (also 
at random) by considering any other of its q admissible values. 

2. This move, from an initial state with energy E and variation 6E, is ac- 
cepted in accordance with the transition probability: 

W{E^ E + AE; Pb) = min {exp [~Pb5E] , 1} . (48) 
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A MC step is produced after considering L'^ moves regardless of whether they 
have been accepted or rejected. 

The extension of this local algorithm with the consideration of an environ- 
ment associated with an arbitrary probabilistic weight oj{E) is achieved by using 
the transition probability Clearly, Eq.® fulfils the detailed balance con- 
dition ([55)) . When the system size A'^ is sufficiently large, such a transition 
probability is practically given by Eq.Q, which is exactly the canonical transi- 
tion probability (j48p modified by the inclusion of a variable inverse temperature, 
/3b — >■ {E). Remarkably, the transition probability ^ can be exactly rewrit- 
ten in the form (|^ by using the value P* = Pui[E,n) corresponding to 6* = 
for the particular case of the linear ansatz (P^ . This is a useful representation 
because both the initial and the final system configurations, X and X\ are a 
priori known for this local MC method. 

3.2.2 Cluster algorithm 

Other important MC methods are the so-called cluster algorithms^ which are 
usually more efficient than any local MC method (Metropolis). However, the 
success of these methods is not universal because the proper cluster moves 
needed seem to be highly dependent on the system, and efficient cluster MC 
methods have only been found for a small number of models [HI [131 [H Ull [23l 

[21[lSl[2i[13- 

The idea of using nonlocal moves was first suggested by Swendsen-Wang 
[m [13] for the case of Ising model and its generalizations, the Potts models. 
Such cluster algorithms are based on a mapping of this model system to a 
random cluster model of percolation throughout the equation: 

Z {Pb) = = J2 (1 - Pf'^' > (49) 

{a} {n} 

where p — p{I3b) = 1 — e~^^, b is the number of bonds and Nc the number 
of clusters. We shall consider in the present study the Swendsen-Wang cluster 
algorithm, whose scheme reads as follows: 

1. Examine each nearest neighbor pair and create a bond with probability 
p (/3b). That is, if the two nearest neighbor spins are the same, a bond is 
created between them with probability p (/3_b); if spin values are different, 
there will be no bond. 

2. Identify clusters as a set of sites connected by zero or more bonds (i.e., 
connected component of a graph) . Relabel each cluster with a fresh new 
value at random. 

The extension of this cluster MC method by using the present ideas is 
achieved by introducing a third step: 

3. Redefine the inverse temperature of the bath fis — Puj [Ei) employed to 
obtain the next system configuration Xi^i by using the energy Ei = E{Xi) 
of the present configuration Xi. 
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Figure 2: Microcanonical caloric and curvature curves, f3 (e) (backbending de- 
pendencies) and k{s) ([/-like dependencies), obtained from MC simulations by 
using the extended version of Metropolis importance sample algorithm for some 
realizations of 2D g-state Potts models. 



While the bath inverse temperature /3b is redefined after every local move 
in the Metropolis importance sample, such a redefinition only takes place in a 
clusters algorithm as the Swendsen-Wang after every MC step because the clus- 
ters moves demand a constancy of the bath temperature. The present method 
is a simple example of extended canonical MC algorithm that does not ful- 
fil the detailed balance condition ([5^ due to its using the approximated in- 
verse temperature = Pu{Ei) instead of the exact value /3* = l3i^{E^) with 
E'Y^ ~ [Ei + ii'i+i] /2. As naturally expected, the violation of the detailed bal- 
ance condition introduces finite size effects in the calculation of the expectation 
values of the thermal fluctuations, which affects the calculation of the heat ca- 
pacity C via Eq. ipTj) . We have verified by mean of preliminary calculations that 
the use of the transition inverse temperature Ptii^ili) instead of the instanta- 
neous value l3u{Ei) in the MC estimation of the expectation values: 

M ^ M 

i=l i=l 

significantly reduces the incidence of such undesirable errors. 
3.3 Results and discussions 

Results of extensive MC calculations by using the extended version of Metropolis 
importance sample (hereafter referred to as extended MIS) are shown in FIGlU 
We are limited to consider here the 2D Potts models with L = 25 for several 
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Figure 3: Comparison between the caloric curves of the 2D 7-state Potts model 
with L = 25 obtained from MC simulations by using the canonical Metropolis 
importance sample (canonical MIS) and its extended version (extended MIS). 
The distribution functions as well as the evolution of the average system energy 
(e) shown in the inset panel correspond to MC runs by using these methods 
with /3e = 1.682. See the text for further explanations. 



values of g, where each point of these curves has been obtained from a MC run 
with 10^ steps. While the microcanonical curvature curve k (e) for the case 
g = 4 only slightly touches the horizontal line with k = 0, the other cases with 
g > 4 clearly exhibit negative values k (e) < 0, that is, macrostates with negative 
heat capacities C < 0. This last observation evidences that the 2D Potts model 
with q — 4: undergoes a continuous phase transition, while the phase transition 
is discontinuous (first-order) for those cases with g > 4 [28] . 

As expected, the anomalous macrostates with C < cannot be accessed by 
using the ordinary Metropolis importance sample (hereafter referred to as canon- 
ical MIS). This limitation is explicitly shown in FIGl3]for the particular case of 
the 2D 7-state Potts model. These results constitute a simple exemplification 
of the schematic behavior represented in FIG[T] For a better understanding, 
we also show here the corresponding energy distribution functions and the dy- 
namical evolutions of the average system energy (e) (inset panel) obtained from 
MC runs by considering both canonical MIS and extended MIS algorithms with 
Pe = 1-682. 

The stationary macrostates at Ea and Ec derived from the intersection of the 
microcanonical curve /3 (e) and the bath inverse temperature /3s — (3^ = 1.682 
are thermodynamically stable within the canonical ensemble. Consequently, the 
canonical energy distribution function is bimodal and its peaks are related to 
the coexisting phases. As consequence of such bimodality, the system energy 
exhibits eventual random transitions between the coexisting peaks, which lead 
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Figure 4: Comparison among the caloric curves of the 2D 7-state Potts model 
with L = 25 obtained from MC simulations using the canonical Swendsen-Wang 
clusters algorithm (canonical SW) and its extended version (extended SW), as 
well as the extended Metropolis importance sample method (extended MIS). 
The energy distribution functions and the dynamical evolution of the average 
system energy (e) shown in the inset panel correspond to MC runs by using the 
above clusters algorithms with Pe — 1-29. 



to a slow equilibration of the corresponding average energy (e) (inset panel) . 

The thermodynamic behavior radically changes when this system is put in a 
thermodynamic situation with non- vanishing correlated fluctuations {6/3ujSE) ^ 
0, which is implemented here by considering a dependence Piu{s) = f5e + X{s—ee), 
where /3e — 1.682, A = Aopt(Kfc) 7^ and Sg = £&, with ki, being the curvature 
at the stationary point located within the anomalous region with C < 0. 
The canonically stable stationary points Sa and Sc become thermodynamically 
unstable and their corresponding peaks disappear from the energy distribution 
function. Conversely, the canonically unstable stationary point Sb now becomes 
thermodynamically stable, and its position determines the maximum of the 
only peak of the energy distribution function. Once the bimodal character of 
the energy distribution function is suppressed, the average system energy (e) 
shows fast convergence towards its equilibrium value £& (inset panel). 

A comparative study of a 2D 7-state Potts model with L = 25 by using the 
Swendsen-Wang cluster algorithm is shown in FIGH) Although the canonical 
Swendsen-Wang method (hereafter referred to as canonical SW) is more efhcient 
than the canonical MIS, it is unable to describe the existence of an anomalous 
region with C < and suffers from a slow relaxation in this kind of situation, 
as it is clearly shown in FIG HI Such limitations are circumvented by using 
its extended version (hereafter referred to as extended SW). As the previously 
discussed example, the extended SW algorithm eliminates the bimodality of the 
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Figure 5: Comparative study between the extended MIS and SW algorithms 
with the known Wang-Landau samphng method. Despite of the relative good 
agreement, one can appreciate some small discrepancies, which could be mainly 
associated with size effects. 



energy distribution function for /3e = 1.29 and leads to a fast convergence of 
the average energy per particle along the simulation. Note also the remarkable 
agreement between the extended MIS and extended SW algorithms despite the 
latter one does not fulfil the detailed balance condition ((39|) . This is a clear 
evidence that a suitable approximation of the inverse temperature of Eq. (|^^ 
is enough to provide a precise estimation of the microcanonical caloric curve 
/3{E) as long as the system size is sufficiently large. 

The comparison among the above extended canonical MC algorithms and 
the known Wang-Landau sampling method flSj is shown in FIGl5] Although 
these results exhibit very good consistency, one can note small but appreciable 
discrepancies within the anomalous region with C < 0. These relative differ- 
ences are naturally expected due to two reasons. Firstly, a direct numerical 
differentiation of the entropy S = log 17 (inset panel) obtained from the Wang- 
Landau method enhances the underlying statistical errors of its MC calculation, 
an therefore, the final result crucially depends on how one defines this math- 
ematical operation for this discrete observable. On the other hand, the work 
equation ([^S]) is supported by the Gaussian character of the energy distribution 
function p{e), which arises as an asymptotic distribution when the system size 
N is sufficiently large. Clearly, the distribution function undergoes small de- 
viations from the Gaussian shape when N is not as large. We shall show in a 
forthcoming paper that these size effects can be taken into account to improve 
the precision of some extended canonical MC algorithms. 

Overall, it is important to remark that the present proposal is focussed on 
the solution of the slow sampling problems observed in large-scale canonical MC 
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Figure 6: Microcanonical caloric curves of the 10-state Potts models for lattice 
dimensions d = (2,3,4) with a fixed number of lattice sites N = L'^ = 4096, 
which were obtained from MC simulations by using the extended SW algorithm 
with 2 X 10^ MC steps for each calculated point. The plot is performed in terms 
of the re-scaled variables e"^ = 2e/d and /3'^ = Pd/2. Inset: Comparison between 
extended SW and Wang-Landau method for the 2D 10-state Potts model with 
L = 64. 



simulations, that is, in systems undergoing a temperature driven discontinuous 
phase transition with sizes sufficiently large to support the Gaussian approxi- 
mation. In particular, the previous extended canonical MC algorithms can be 
useful to study Potts models with many spin states and higher dimensions, as 
the cases shown in FIGlHl where we illustrate the caloric curves derived from 
MC simulations by using the extended SW algorithm for q — 10^ d = (2,3,4), 
and a fixed number of lattice sites N = L'^ = 4096. The comparative study case 
with the Wang-Landau method shown in the inset panel allows us to verify that 
the agreement between these methods is more significant with the increasing of 
N. 

To quantitatively characterize the efficiency of the previously discussed ex- 
tended canonical MC methods, one should obtain the decorrelation time t, that 
is, the minimum number of MC steps needed to generate effectively indepen- 
dent, identically distributed samples in the Markov chain. Its calculation is 
performed here by using the expression: 

r r M-var{eM) 
T — hm tm — lim ; — r — , (51) 

M^oo M^oo var (ei) 

where var (em) = (^If ) — (em)^ is the variance of Em, which is defined as the 
arithmetic mean of the energy per particle e over M samples (consecutive MC 
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steps): 

1 

^M = ^E^- (52) 

i=l 

This quantity is calculated for the particular MC runs considered in FIG [3] and 
FIGIH This study allows us to verify that the use of a variable dependence 
I3^{e) instead of a constant parameter Pb enables the reduction of the decor- 
relation time of the Metropolis importance sample from r ~ 450 to r ~ 80. 
The improvement is even more significant for the Swendsen-Wang clusters al- 
gorithm, which experiences a reduction of the decorrelation time from r ~ 300 
to r ~ 12.2. 

Let us now analyze the behavior of the decorrelation time t with the increas- 
ing system size N = L"^ at the inverse temperature (3c of the first-order phase 
transition. Such a study demands the performance of preliminary calculations 
of the quantity /3j, due to its underlying dependence on the system size N. For 
computational limitations, we decide to restrict our analysis for the case of 2D 
10-state Potts model with L = (8, 16, 32, 64), whose results are shown in FIGl?! 
The entropy per particle s (e) was obtained from the second-order approximation 
of the power-expansion As = (3Ae— ^^Ae^, which allows us a direct calculation 
of this thermodynamic function through the inverse temperature /3 and curva- 
ture K. These numerical results are shown in panel b) of FIG (71 or more exactly, 
the quantity s* (e) = s {e) — j3ce — sq, with sq being a suitable constant, which 
allows us to appreciate better the convex intruder of the entropy accounting 
for the existence of macrostates with C < 0. The mathematical form of this 
last anomaly enables the acquisition of the latent heat ql — £3 ~ £1 and the 
entropy-loss per particle Assurf associated with the existence of surface correla- 
tions [3], and (£1,62, £3) represent the three stationary points where the caloric 
curve /3 (£) takes the value of the phase transition inverse temperature /3c- It 
worth to clarify that the critical value /3c is simply the point of discontinuity of 
the first derivative of the Planck thermodynamic potential p(/3) estimated from 
the microcanonical entropy s{e) via the Legrendre's transformation; 

p(/3) = min[/3£-s(e)]. (53) 

In fact, once obtained the microcanonical entropy S{E), one can calculate any 
thermo-statistical quantity by using the canonical distribution function: 

dp{E\l3) = Ae-^'^+^^'^'^dE. (54) 

Relevant physical observables and thermodynamic parameters derived from this 
type of analysis are summarized in Table [TJ 

Results of extensive calculations of dependencies of the decorrelation time r 
versus the system size A'' at the point of phase transition /3c are shown in F1G|51 
It is clearly evident that the extended MC methods are always much more effi- 
cient than their respective canonical counterparts. In particular, the extended 
SW method reduces the exponential growth of the decorrelation time r with TV 
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L 


8 


16 


32 


64 




1.415 


1.422 


1.424 


1.426 


£i 


0.319 


0.319 


0.321 


0.329 


£2 


0.767 


0.755 


0.737 


0.719 


£3 


1.165 


1.114 


1.074 


1.049 


QL 


0.846 


0.795 


0.753 


0.72 


ASsurf X 10^ 


72 


5.6 


3.1 


1.7 



Table 1: Dependence of some physical observables and thermodynamic param- 
eters of the 2D 10-states Potts model with the increasing lattice size L. 

to a weak power-law dependence r oc with a = 0.2. As expected, the ex- 
tended MIS is less efficient than the extended SW. However, the iV-dependence 
of its corresponding decorrelation time r does not differ in a significant way 
from the one associated with the cluster algorithm. Since the calculation of 
the decorrelation time r of the canonical MC methods demands very large MC 
runs, we use the 1/M extrapolation in order to obtain some rough estimates of 
the non-equilibrated points. 

According to Berg [29j, the multicanonical method and its variant are able 
to reduce the exponential divergence of the decorrelation time r to a power 
dependence with typical exponent a = 2 ~ 2.5. Although such a power-law 
behavior accounts for a less efficient convergence than the one achieved by using 
present proposal, it is important to remark that the multicanonical methods 
possibilities an effective exploration of the entire energy region in a single MC 
run. Conversely, our methodology only explores a small energy region (se — 
a, Ee + cr) with a typical width a ~ so that, a minimum of n oc MC 

runs are required to study the same energy range of multicanonical methods. 
Thus, the size dependence the number of MC samples M can be estimated by 
considering the number of MC runs n and the decorrelation time r as follows: 

M~nr-7V"', (55) 

which leads to the effective exponent a* ~ 0.7. Such an effective power-law 
growth still considers a more efficient convergence rate than the one associated 
with multicanonical methods, overall, when one also takes into account the fact 
that these re-weighting techniques demand a preliminary reconstruction of the 
probabilistic weight i-o{E), which is a procedure that consumes a significant 
fraction of the available computational resources. 

Let us finally reconsider the question about the optimal value of the cou- 
pling constant A that allows to achieve the best efficiency by using an extended 
canonical MC algorithm. An example of such a study is shown in FIG El where 
we illustrate the behavior of the decorrelation time r and the total dispersion 
A|, — {AEY / N -\- N {APi^Y for positive values of the coupling constant A. The 
data was obtained from MC simulations of the 2D 10-state Potts model with 
L = 32 by using the extended SW clusters algorithm with environment inverse 
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temperature Pu:{e) — j3e + A(e — £&), whose parameters /3e = j3c and Eg = £2 in 
order to access the canonically unstable stationary macrostate 82 located within 
the anomalous region with C < at the point of the first-order phase transition 

As expected, the decorrelation time r and the total dispersion show 
large values when the coupling constant A is close to zero as a consequence of 
the imposition of the external conditions associated with the canonical ensemble 
(m , where take place thermodynamic instability of the macrostates with C < 
and the bimodal character of the energy distribution function p{e). The increase 
of the coupling parameter A for values that obeys the condition ([55]) produces a 
progressive reduction of the energy dispersion /S.E in Eq. ([5^ . This last behavior 
leads to a reduction of the decorrelation time r due to the MC sampling performs 
a faster exploration for a smaller energetic region. However, the decorrelation 
time T starts to grow after reach its minimum value at a certain At since the 
system configurations X cannot be modified in an appreciable way by the MC 
sampling if the energy is only allowed to undergo too small changes around its 
equilibrium value. 

For a general case of the environment inverse temperature (|29)) with A 7^ 0, 
both the system fluctuating behavior and the decorrelation time depend on the 
particular value of the coupling constant A. In such cases, the efficiency of 
an extended canonical MC method can be evaluated throughout the minimal 
number of MC steps needed to achieve a certain precision in the calculation of a 
given thermo-statistical observable. The statistical error Se associated with a set 
of n independent outcomes {xi, i = 1, ■■n} can be estimated from the standard 
deviation a as 5e^ = jn. In MC calculations, independent samples are only 
obtained after considering a number of MC steps equal to the decorrelation 
time T. Thus, the number of MC steps S needed to obtain an estimation of a 
given point (e, /?) of the caloric curve with a precision \J Se^ + (5/3^ < 8a can be 
evaluated in terms of the decorrelation time r, the total dispersion and the 
system size TV as follows: 

* " - ^^^^ 

The quantity 77 = tA|, could be referred to as the efficiency factor. Clearly, a 
dynamic criterion in order to provide an optimal value of the coupling constant 
A to achieve the best efficiency by using this kind of MC calculations is given 
by minimizing the efficiency factor 77. 

According to the data shown in FIGEl the value of the coupling constant 
A^ corresponding to the minimum of the efficiency factor 77 is located within 
the interval between the position of the minimum At of the decorrelation time 
r and the minimum As of the total dispersion A^, A € [As,At]. Remarkably, 
the efficiency factor 77 does not change in a significant way within this last 
interval. By considering this last observation and the fact that the calculation 
of the correlation time r demands extensive calculations, one can realize that 
the value of the coupling constant As corresponding to the minimum of the total 
dispersion A|,, Eq . ([BS]) . provides a good value in order to perform a very efficient 
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MC calculation of the caloric curve by using extended canonical algorithms. 
In general, it is always convenient to keep as small as possible the value of 
the coupling constant A. Indeed, the statistical error e involved in the MC 
calculation of the standard deviation a1 = (AE)'^ /N considered to obtain the 
system curvature k — f3'^N/C from Eq. (PT|) leads to the existence of a statistical 
error 6k ^ (A+1)(A+K)e, which grows with the increase of the coupling constant 
A. 

4 Final remarks 

We have shown that conventional Monte Carlo methods based on the consid- 
eration of the Gibbs canonical ensemble (jlj can be easily extended in order to 
capture the existence of an anomalous regime with negative heat capacities and 
avoid the incidence of the super-critical slowing down. The key ingredient is 
to replace the use of a bath with an infinite heat capacity by an environment 
with a finite heat capacity that obeys Thirring's constrain (|22p . Such an equi- 
librium situation, characterized by the existence of non-vanishing correlations 
{5I3^5E) ^ 0, is inspired on the generalized equilibrium fluctuation-dissipation 
relation ([2]), which allows to introduce several improvements to the methodol- 
ogy of Gerling and Hiiller based on the consideration of the dynamical ensemble 
(ESI). 

The way to introduce an environment with a variable inverse temperature 
depends on the own features of each canonical Monte Carlo method, although 
such a question seems not to be a difficult problem in the case of classical al- 
gorithms. While it could be desirable that the implementation of this kind of 
methodology obeyed the detailed balance condition ([5^ . the application ex- 
amples considered in the section [3] show that one can still obtain a good MC 
estimation of the microcanonical caloric curve I3{E) without fulfilling the de- 
tailed balance as long as the system under analysis be sufficiently large. 

Before concluding this section, it is worthwhile to mention that Eq.(l2]) con- 
stitutes a particular case of a more general equilibrium fluctuation-dissipation 
theorem, which accounts for the system fluctuating behavior in a thermody- 
namic situation characterized by the incidence of several control parameters 
[30] . Roughly speaking, this theorem provides a general extension of some 
other well-known fluctuation relations such as the one involving the isothermal 
magnetic susceptibility and magnetization fluctuations of a magnetic system, 
XT — f3{SM'^'^, or the isothermal compressibility and volume fluctuations of 
a fluid system, VKt — /3((5V^^), which are compatible with the existence of 
anomalous response functions, e.g., negative isothermal susceptibilities xt < 
or negative isothermal compressibilities Kt < 0. Clearly, this general theorem 
suggests a direct extension of the present methodology in order to enhance MC 
methods based on the so-called Boltzmann- Gibbs distributions: 

dpBG{E, X) = ^ exp [~13{E + YX)] dEdX (57) 
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to account for the existence of macrostates with anomalous response functions. 
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Figure 7: Study of the thermodynamic behavior of the 2D 10-state Potts model 
for different lattice sizes L: Panel a) microcanonical caloric curves /3 (e), Panel 
b) microcanonical entropy plotted in term of the quantity s* (e) = s (e) — ^cE— So 
in order to reveal the convex intruder associated with the existence of surface 
correlations. 
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Figure 8: Dependencies of the dccorrclation time r versus the system size N at 
tlie point of first-order phase transition /3c for different MC methods considered 
in the study of the 2D 10-state Potts model. Inset panel: A more detailed plot 
for the decorrelation time of the extended SW method, which clearly shows a 
weak power-law dependence r oc N"' with a = 0.2. 
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Figure 9: Dependence of the decorrelation time r, the total dispersion and 
the efficiency factor t] = tA|, on the coupling constant A obtained from MC 
simulations of the 2D 10-state Potts model with i = 32 by using the extended 
SW clusters algorithm. See the text for further explanations. 
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